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Abstract 

In order to generate initial data for nonlinear relativistic simulations, one needs 
to solve the Einstein constraints, which can be cast into a coupled set of nonlin- 
ear elliptic equations. Here we present an approach for solving these equations on 
three-dimensional multi-block domains using finite element methods. We illustrate 
our approach on a simple example of Brill wave initial data, with the constraints 
reducing to a single linear elliptic equation for the conformal factor ip. We use 
quadratic Lagrange elements on semi-structured simplicial meshes, obtained by tri- 
angulation of multi-block grids. In the case of uniform refinement the scheme is 
super convergent at most mesh vertices, due to local symmetry of the finite element 
basis with respect to local spatial inversions. We show that in the superconvergent 
case subsequent unstructured mesh refinements do not improve the quality of our 
initial data. As proof of concept that this approach is feasible for generating multi- 
block initial data in three dimensions, after constructing the initial data we evolve 
them in time using a high order finite-differencing multi-block approach and extract 
the gravitational waves from the numerical solution. 



1 Introduction 

This paper is part of an effort to numerically solve Einstein's equations on general domains 
with arbitrary shapes and topologies. Accurate simulations on such domains require multiple 
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numerical grids covering different areas of the domain, possibly refined where the solution 
has interesting features. For the numerical discretization one can use either finite differences 
(FD), spectral methods, or finite elements (FE). Current research in numerical relativity 
community is dominated by finite differences methods, not the last reason for that being 
the relative ease of finite differences implementation and parallelization. A lot of important 
results in areas of hyperbolic solvers, boundary conditions, domain representation, appar- 
ent horizon finders, wave extraction, stability proofs etc., were established and successfully 
implemented using the FD methods (see pl2] and references therein). 

Our approach was designed to allow for fast and accurate transformation of the numerical 
function between FE and FD representations. We demonstrate that the FE solution to an 
elliptic equation, generated for the given FD grid, is sufficiently accurate and convergent for 
3D relativistic simulations with independent finite difference hyperbolic solver (QUILT). 

Spectral methods are also widely used [3lH] . mainly for solving elliptic equations, such as 
those arising in initial value problems [SIEEISE], constrained evolution [TU] and apparent 
horizon finders [TT][T^ . They can as well be applied to hyperbolic problems [13], as in the 
recent very successful simulations of the binary black holes merger [13] (using an approach, 
developed by the Cornell- Caltech collaboration [T5lll6lll7|[T8] ). Spectral methods can produce 
exceptionally accurate results, but in order to use these results in a FD evolution code 
one needs to perform an additional interpolation step. The interpolation can be too slow 
for some applications, for instance, in constrained evolution schemes. In our approach, the 
interpolation is trivial, and the number of degrees of freedom on FE and FD grid is the 
same, which allows to exchange solution between the FE and FD solvers without reducing its 
accuracy by interpolation. Another advantage of FE methods compared to spectral methods 
is that the matrices resulting from FE discretizations are sparse. 

Finite element methods can also be applied to hyperbolic problems, however, stable and 
accurate discretizations require using discontinuous Galerkin (DC) flavor of finite elements 
(otherwise the approximation is suboptimal: order of convergence falls to {p — 1) and the 
solution is not necessarily stable). 




Fig. 1. Types of FE meshes. Top left: unstructured mesh, top right: regular mesh, bottom: 
semistructured meshes. 

In this paper, the term grid is used in reference to the finite differences grid (an ordered 
set of points), as opposed to the term mesh for the finite element triangulation. Domain tri- 
angulations can be classified as regular, unstructured or semistructured, depending on their 
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degree of symmetry (see figured]). Finite element methods can handle domains of arbitrarily 
complex shapes and topologies, since every geometrical shape admits approximation with 
succesively refined unstructured symplectic triangulations. However, the domains encoun- 
tered in numerical relativity are relatively simple in the sense that they can be covered with 
semistructured meshes. We argue in favor of using semistructured meshes over unstructured 
ones. One of the motivations is that they can be quickly constructed on top of existing 
FD grids without the need to explicitly maintain all the information about mesh elements, 
such as edges, faces etc. Raising the order of finite elements, one can essentially achieve the 
spectral convergence; however, here we concentrate on lower-order finite elements (first- and 
second-order), and show how the superconvergence on semistructured grids for second-order 
elements can be exploited to obtain accurate 4-th order convergent solution, suitable for use 
in finite differences simulation. 

The method is tested on the particular example of multipatch FD grid. In differential ge- 
ometry it is common to define a manifold through a set of possibly overlapping patches, 
with each patch mapped into an open, simply-connected subset of Euclidean space (see, for 
example, [IS]). This is a natural way of describing manifolds with nontrivial topology, which 
cannot be covered by a single coordinate chart. After discretization every chart becomes a 
FD grid, which can be thought of as a discrete coordinate chart, mapped into a region of . 
At the continuum level, the patches are glued together by coordinate transformations in the 
areas where they intersect. At the discrete level, if the neighboring patches do not have com- 
mon points interpolation can be used to fill in the missing ones; this approach is commonly 
referred to as overlapping- grids. On the other hand, in a multi-block approach the patches 
abut rather than overlap, and the grids are constructed in such a way that neighboring grids 
share boundary points. 

A multipatch approach in numerical relativity has several advantages. In many cases the do- 
main of interest is asymptotically fiat. If it contains one or more black holes, then each hole 
can be excised from the computational domain by introducing an inner smooth boundary 
around the singularity. If appropriately chosen, this boundary does not require any physical 
boundary conditions, since all characteristic modes leave the domain. It is preferable that 
this boundary is smooth, which in general requires the use of multiple patches. Similarly, the 
preferred shape for the outer boundary when simulating asymptotically fiat spacetimes is a 
sphere, as this is the topology of null infinity J'~^, which is best suited for extracting gravi- 
tational waves. A multipatch domain structure easily accomodates both type of boundaries 
while avoiding coordinate singularities, such as those associated with spherical or cylindrical 
coordinates. The use of multiple patches is unavoidable in cosmological simulations with 
non-trivial topologies. In addition, multiple patches are in general more efficient than regu- 
lar grids, since they decouple different spatial directions. For example, under conditions of 
approximate spherical symmetry, one could surround the system with spherical patches and 
fix the number of points in angular direction, while increasing the domain only in radial 
direction. This way, pushing the outer boundary out becomes an order problem, as op- 
posed to A^^ with Cartesian grids. This same feature makes them useful, in particular, for 
many relativistic astrophysical studies which are assumed to be approximately spherically 
symmetric [20] . 
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Einstein's equations are often written in a form such that the equations divide into hyperbohc 
(evolution) and elhptic (constraint) sub-systems. When solving the Einstein vacuum (that 
is, in the absence of matter) evolution equations, these can be cast in linearly-degenerate 
form, and the solutions are in general expected to be smooth. Those cases are ideally suited 
for high order or spectral methods. When including matter, on the other hand, for example 
when dealing with the general relativistic hydrodynamics equations, shocks are expected. In 
those cases a possible choice is high resolution shock capturing methods [2T|22] adapted to 
the presence of several patches for the hydrodynamical sub-system and high order or spectral 
methods for the metric sub-system [20123] • 

This paper is concerned with the elliptic sector of Einstein's equations, more precisely with 
the generation of initial data needed for evolutions on multipatch geometries. Initial data 
on a spacelike hyperslice has to satisfy a set of Enistein's constraints, which can be cast 
into a coupled nonlinear elliptic system of partial differential equations (PDE). This el- 
liptic PDE system has been studied extensively for many years, with a complete solution 
theory developed in the case of domains which are closed (i.e. compact without bound- 
ary) 3-manifold spatial slices of spacetime with constant or nearly-constant mean extrinsic 
curvature p^f25] . More recently, this solution theory has been extended to both closed 3- 
manifolds and to compact 3-manifolds with boundary, having mean extrinsic curvature far 
from constant [2B|27f28j . It has also been shown recently [2^3UII31j that geometric partial 
differential equations of this type can be solved accurately and efficiently using adaptive fi- 
nite element methods on manifolds with general topologies. For general multi-block systems, 
the additional structure allows one to construct semi-structured triangulations, which gen- 
erate super convergent finite element solutions. The property of superconvergence holds at 
the vertices of the multi-block grid triangulation, which simplifies transport of the solution 
to the finite difference grid, since no interpolation is required. 

The above considerations define the approach of this paper, which involves solving the Ein- 
stein constraint equations on semi-structured multi-block triangulations using finite element 
methods. We apply this approach to the Brill wave initial data problem, which represents 
a linear vacuum case with zero extrinsic curvature, with the constraints system reducing 
to a single elliptic equation Aip — Vip = for the conformal factor ip (where A is the fiat 
3-dimensional laplacian and V is a potential function related to the Ricci scalar of the freely 
given metric). Nevertheless, the results presented are expected to hold in general case, be- 
cause the property of superconvergence is the property of the mesh and the finite element 
spaces used. If this property holds for a hnearized problem, it is also expected to hold for 
the full nonlinear one (see Chapter 9 of |32j). 

The paper is organized as follows. Section [2] gives an overview and summary of the Finite 
Element Toolkit (FETK) [29], which we here use for solving the Einstein constraint equations 
with finite element methods. Section [3] discusses our approach to multi-block triangulations 
and, in particular, why superconvergence is expected. In section H] we evaluate the accuracy 
of our approach by solving for several three-dimensional elliptic test problems on multi-block 
domains with known exact solutions. Finally, in section \5\ we solve the Einstein constraints 
on a multi-block domain for the case of Brill waves. After constructing initial data, we use 
the multi-patch infrastructure QUILT to solve for the Einstein evolution equations in time 
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and extract gravitational waves from the numerical solution. QUILT is an ongoing effort, 
more details of which can be found in [33], [31], [35], [36], [37], [20]. In tensor notations, used 
throughout the paper, Latin indices from the beginning of the alphabet (a, 6,...) run from 
to 3, Latin indices i, j, etc. run from 1 to 3, and the usual Einstein summation rule is 
assumed. 



2 Solving elliptic PDEs using the Finite Element ToolKit 

In this paper we use the Finite Element ToolKit (FETK) [29] (see also [38)131] ) to solve the 
Einstein constraint equations. FETK is an adaptive multilevel finite element code developed 
over a number of years by the Hoist research group at UC San Diego and their collaborators 
(see also [39]). It is designed to produce provably accurate numerical solutions to nonlin- 
ear elliptic systems of tensor equations on (Riemannian) multi-dimensional manifolds in an 
optimal or nearly-optimal way. We will summarize the main features of FETK here; more 
detailed discussions of its use for general geometric PDF may be found in [29|30] . and specific 
application to the Einstein constraint equations may be found in [3T]|^ . 

FETK contains an implementation of a "solve-estimate-refine" algorithm, employing inexact 
Newton iterations to treat non-linearities. The linear Newton equations at each inexact 
Newton iteration are solved with unstructured algebraic multilevel methods which have been 
constructured to have optimal or near-optimal space and time complexity (see [11112]). The 
algorithm is supplemented with a continuation technique when necessary. FETK employs a 
posteriori error estimation and adaptive simplex subdivision to produce provably convergent 
adaptive solutions (see [13111]). 

Several of the features of FETK are somewhat unusual, some of which are: 

• Abstraction of the elliptic system: The elliptic system is defined only through a nonlinear 
weak form along with an associated linearization form over the domain manifold. To use 
the a posteriori error estimator, a third function F{u) must also be provided (essentially 
the strong form of the problem). 

• Abstraction of the domain manifold: The domain manifold is specified by giving a polyhe- 
dral representation of the topology, along with an abstract set of coordinate labels of the 
user's interpretation, possibly consisting of multiple charts. FETK works only with the 
topology of the domain, the connectivity of the polyhedral representation. 

• Dimension independence: The same code paths are taken for two-, three- and higher- 
dimensional problems. To achieve this dimension independence, FETK employs the sim- 
plex as its fundamental geometrical object for defining finite element bases. 

2. 1 Weak Formulation Example 

We give a simple example to illustrate how to construct a weak formulation of a given PDF. 
Here we assume the 3-metric to be flat so that V is the ordinary gradient operator and ■ the 
usual inner product. Let AA represent a connected domain in M? with a smooth orientable 
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boundary dAi, formed from two disjoint 2-dimensional surfaces doAi and diM. 
Consider now the following semilinear equation on Ai: 



~V ■ {a{x)Vu{x)) +b{x,u{x))=0 in M, (2.1) 

n{x) ■ {a{x)Vu{x)) +c{x,u{x)) = on diM, (2.2) 

u(x) = f{x) ondoM, (2.3) 

where n{x) : dAi ^ is the unit normal to dAi, and where 

a-.M^R^''^, b:MxR^R, (2.4) 

c : diM xR^R, f : doM ^ R. (2.5) 



If the boundary function / is regular enough so that / G H^^'^{doAi), then from the Trace 
Theorem [15], there exists u G H^{Ai) such that / = u\qqM in the sense of the Trace op- 
erator (where H^{At) is the Hilbert space of all real L2-integrable functions on At with 
L2-integrable weak derivative [IH].) Employing such a function u G H^{Ai), the weak for- 
mulation has the form: 

Find ueu + Hi o{Ai) s.t. (F(n), v) = 0, ^ vE HIjj{M), (2.6) 

where the nonlinear form is defined as: 

{F{u),v) = / {aWu ■ Wv + b{x,u)v) dx + c{x,u)v ds. (2.7) 

Jm JdiM 

The "weak" formulation of the problem given by equation (12. 6p imposes only one order 
of differentiability on the solution u, and only in the weak or distributional sense. Under 
suitable growth conditions on the nonlinearities b and c, it can be shown that this weak 
formulation makes sense, in that the form {F{-),-) is finite for all arguments, and further 
that there exists a (potentially unique) solution to (12. 6p . In the specific case of the individual 
and coupled Einstein constraint equations, such weak formulations are derived and analyzed 
in fimiT] . 

To analyze linearization stability, or to apply a numerical algorithm such as Newton's 
method, we will need the bilinear linearization form {DF[u)w,v), produced as the formal 
Gateaux derivative of the nonlinear form {F{u),v): 

{DF{u)w,v) = -^{F{u + ew),v) 
= I [aVw ■ Vv + ^^^^''^^ wv] dx+ f ^^!^^wv ds. (2.8) 

JM \ OU J JdiM OU 

Now that the nonlinear weak form {F{u),v) and the associated bilinear linearization form 
{DF{u)w,v) are defined as integrals, they can be evaluated using numerical quadrature to 
assemble a Galerkin-type discretization involving expansion of m in a finite-dimensional basis. 



6 



As was the case for the nonhnear residual (F(-), ■), the matrix representing the bihnear form 
in the Newton iteration is easily assembled, regardless of the complexity of the bilinear form 
{DF{-)-, ■). In particular, the algebraic system for w = J2j=il3j(pj has the form: 

AU = F, Ui = A, (2.9) 

where 

n 

= {DF{uh + J2 (^k<Pk)<t>v ^i)^ (2-10) 

k=l 

n 

F, = {F{uh + Y.»^)^^^)■ (2-11) 

and {(f)i}fLi, {ipj}jLi are the bases of the test and trial iV-dimensional finite element spaces. 
As long as the integral forms (F(-), ■) and {DF{-)-, ■) can be evaluated at individual points 
in the domain, then quadrature can be used to build the Newton equations, regardless of 
the complexity of the forms. This is one of the most powerful features of the finite element 
method, and is exploited by FETK to make possible the representation and discretization of 
very general geometric PDEs on manifolds. It should be noted that there is a subtle difference 
between the approach outlined here (typical for a nonlinear finite element approximation) 
and that usually taken when applying a Newton-iteration to a nonlinear finite difference ap- 
proximation. In particular, in the finite difference setting the discrete equations are linearized 
explicitly by computing the Jacobian of the system of nonlinear algebraic equations. In the 
finite element setting, the commutativity of linearization and discretization is exploited; the 
Newton iteration is actually performed in function space, with discretization occurring "at 
the last moment." 



3 Generating semi-structured multi-block triangulations with superconvergent 
properties 

In order to avoid the extra step of interpolating the finite element numerical solution to a 
multi-block grid, we use finite element meshes with vertices located at the multi-block grid 
points. We construct such meshes by dividing a convex hull of the set of grid points into 
simplices with vertices only at those points. This procedure of building a simplicial mesh 
based on a set of points is usually referred to as triangulation. Delaunay's triangulation (see, 
for instance, [17]) is an example of such a procedure, generating meshes with simplices of 
the highest possible quality in flat space. Although Delaunay's triangulation minimizes a 
condition number in a resulting linear system, it is mostly used for sets of points of general 
position since it generates completely unstructured meshes; for more regular grids its algo- 
rithm is unnecessarily complex. Here we use a simpler and more straightforward algorithm 
to generate semi-structured meshes with the additional advantage of having superconvergent 
properties. 

The term superconvergence applies in various contexts, when the local or global convergence 
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order of a numerical solution is higher than one would expect P8|49|32|[50] . The type of 
superconvergence we are interested in is superconvergence by local symmetry [32j, which 
occurs in function values for discretizations of second order elliptic boundary value problems, 
and amounts to an additional < a < 1 in the convergence order. It occurs at any given 
point in which the finite element basis is locally symmetric, or approximately symmetric, 
with respect to local spatial inversion at that point (see figure [2]). The triangulation method 
that we use produces simplicial meshes with this type of symmetry at all vertices inside 
each block, and at some vertices at the interblock boundaries. Therefore, in our meshes we 
expect superconvergence everywhere but at the non-symmetric interblock boundary points. 
For second order elliptic equations and piecewise polynomial finite elements of even degree 
superconvergence occurs in the solution itself, while for odd degree it occurs in the first 
derivatives of the solution [5T][52]. 



Fig. 2. Example of two-dimensional simplectic mesh, locally symmetric with respect to the point 

m- 

With this in view, the quadratic Lagrange finite elements at the nodal points of the meshes 
that we use are expected to give solutions with convergence order 3 + cr (for some < tr < 1). 
This is an important advantage of multi-block triangulations: the resulting grid solutions have 
a convergence rate which is even higher than the global convergence rate of the original finite 
element solution, almost everywhere. 

We start by dividing a cell, then combine cells into a block, and finally combine blocks to- 
gether in a conforming way. Next we explain how these procedures are performed. The main 
complexity in building multi-block triangulations comes from the last step, since triangula- 
tions of each block are not necessarily conforming at their interfaces. We will show, however, 
that our method guarantees conforming simplicial triangulations. 

Cell triangulation. For a cubical cell, there exist two possible triangulations (see figure [3]): 
the cell can be divided into either five or six tetrahedra. The former is known as middle cut 
triangulation [53], and the latter as Kuhn's triangulation [SI]. While Kuhn's triangulation 
produces simplices of equal shape and can be more easily extended to an arbitrary number 
of dimensions [SS], the middle cut triangulation produces higher quality tetrahedra, in the 
sense that their angles are less acute than in a Kuhn triangulation, which produces to a finite 
element matrix system with better condition number. Usually, the quality of tetrahedra is 
described by the ratio R/3r of radii of circumscribed (i?) to inscribed (r) spheres. For a 
middle cut triangulation the central tetrahedron has, obviously, minimum possible aspect 
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ratio -R/3r = 1, while for corner tetrahedra we have R/?>r = ~ 1-37. At the same time, 
in Kuhn's triangulation all tetrahedra have aspect ratio 1 + l/^/S ~ 1.58. 




Fig. 3. Middle cut triangulation and Kuhn's triangulation: dividing a cube into five or six tetra- 
hedra, respectively. 

Block triangulation. A block can be triangulated in many possible ways. We will consider 
the two most straightforward and symmetric ones, referring to them as uniform and clustered 
block triangulations. Both types of triangulation are locally symmetric with respect to local 
space inversion at any inner vertex of the mesh. In uniform block triangulation, exactly 
the same cell triangulation, with the same orientation, is applied to all cells in a block (see 
figure H]). This produces a conforming mesh for Kuhn's triangulation, but fails to produce a 
conforming one for the middle cut one, since the triangulation patterns between neighboring 
cell interfaces are not compatible. 




Fig. 4. Uniform (left) and clustered (right) block triangulations. With Kuhn's triangulation of each 
cell, both types of block triangulation can be done in conforming way for neighboring cells. Cells, 
triangulated by middle cut, can only be arranged into clustered configuration. 




Fig. 5. Union jack triangulation of 2 x 2 x 2-cell cluster. The triangulation pattern on each side 
of the cluster is the same. Inside, each cluster can have either middle cut (left) or Kuhn (right) 
triangulations. 

To handle this compatibility problem, it is convenient to group neighboring cells into 2x2x2 
clusters and triangulate each cluster as shown in figure O with "union jack" patterns on 
each side. This guarantees conformity between neighboring clusters. We will label this type 
of triangulation as clustered block triangulation. It can be used with both types of cell 
triangulation. 
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Multi-Block triangulation. In a multi-block system, both uniform and clustered triangu- 
lations of each block can be arranged into a conforming simplicial mesh. Obviously, clustered 
triangulations of each block with even number of cells on the interfaces between the blocks 
assemble themselves in a conforming manner. A uniform triangulations arrangement can 
be constructed from clustered triangulation of the same system by first removing and then 
replicating layers of cells with Kuhn's triangulation. Two different types of semi-structured 
multi-block triangulations, used in this paper, are shown in figure [6l 




Fig. 6. Seven-block system for the sphere (left), and its triangulations, generated using clustered 
(center) and uniform (right) block triangulations. 



4 Quality of our finite element solutions on semi-structured multi-block trian- 
gulations 

In this section we perform a numerical study of the quality of our finite element solutions ob- 
tained using semi-structured multi-block triangulations. We investigate not only the solution 
itself, but also answer the question of whether this solution is appropriate for finite-difference 
evolution codes with high-order numerical derivative operators. After introducing the domain 
structure and weak formulation of the second-order elliptic equation, we evaluate the solu- 
tion convergence order, and show superconvergence for quadratic finite elements. Then we 
evaluate the finite element solution at the multi-block gridpoints, apply various high-order 
finite difference operators and check the convergence orders of its first and second numerical 
derivatives. The observed convergence orders are consistent with the expected values. 

4.1 Domain structures 

Both FETK and QUILT were developed to handle equations on general manifolds, with an 
arbitrary number of charts. However, for the domain structures considered in this paper we 
can embed the computational domain into a reference Euclidean 3-dimensional space with a 
single fixed Cartesian system of coordinates to label all vertices and nodes of the mesh, and 
we do so. The domain of interest here will then be a spherical domain of radius R, equipped 
with a seven-blocks or thirteen-blocks system (see figure [7]), with local patch coordinate 
transformations defined as in [55], The seven-block geometry (see figure [7]) is fuHy 
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(a) 



Fig. 7. An equatorial cut of (a) seven-block and (b) thirteen-block systems, (c) Grid dimensions 
for the seven-block system. 

specified by fixing the outer sphere radius Rout and the side of the inner cubical patch ac- In 
the innermost, cubical patch, there are N^xNyXN^ points, where = Ny = := A^. Since 
the grids are conforming, in the six blocks surrounding the inner one there are N x N x Nj. 
points. The thirteen-block geometry, in turn, can be seen as a seven-block one of radius 
Rmed and N X N X Nr^in points, surrounded by six additional blocks with N x N x Nr^out 
points each, with a total radius Rout- The advantage of the thirteen-block setup is that the 
transversal grid layers in the outermost system are perfectly spherical surfaces. These are 
very convenient in certain applications which require integration over such surfaces, as in the 
the study of the multipole structure of radiation in the wave zone (see, for example, [37]). 
since no interpolation is needed for those integrations (resulting in both higher accuracy and 
speed). Throughout this paper, when we change the resolution we keep the ratios : A',, 
and : N^^in : N^^out in the seven- and thirteen-block systems, respectively, fixed. Multiple 
domains with different values of A^ produce a sequence of domain triangulations {T^} with 
maximal simplex diameters h inversely proportional to A^. Therefore, it is convenient to use 
A^ as a scaling factor in convergence tests, and we do so. 



4-2 Second-order elliptic equation and its weak form 

Let Sr represent our spherical domain with radius i?, centered at the origin, and let BSr be 
its outer boundary. The equations of interest in this paper are of the form: 



where ^/'d is a Dirichlet boundary condition, and is the Laplace operator. We only consider 
the case when both the potential V and the boundary conditions ipD are axisymmetric. 
Following the weak formulation example above (section 12.11) . we obtain the nonlinear weak 
form and the bilinear linearization form: 
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(F(V^),0)=/ {Vij-V(j) + ijV(j)) dx (4.3) 
{DF{ij)x, 4>)= f (Vx ■ V0 + xV^) dx (4.4) 

JSr 

with ip{x) E ijj + Hq{Sr), and 0(a;) , xi^) ^ Hq{Sfi) as discussed in section [21 These two 
forms and the Dirichlet boundary function are everything we need to specify the problem 
in FETK. 

We also consider the same problem with more complex Robin boundary conditions: 

d,^{x) = ^—y^ on OSr , (4.5) 
which has the following nonlinear weak and bilinear linearization forms, 



{F{^),(f))= f (ytp -Vcp + ^pV^) dx+ ^ f {i>-l)(j)ds (4.6) 

J Sr 11 J OSr 

{DF{^)X, 0) = / (Vx ■ V0 + xV(p) dx+^l x<pds (4.7) 

JSr K JdSR 



where now '?/;(x), 0(x), x{j^) ^ //^(S'j^). We will use these Robin boundary conditions when 
solving for Brill waves below in section [5l 

4-3 Testing convergence of the solution 

As a test problem for our approach to solving elliptic equations on a semi-structured grid 
using finite elements, we solve equation (14.11) with three different potentials: 

Va = -3a^^ 



2\2 ' 

oJ 



where R = r'^ — cTrv'^ cosh Z = + a1, C = cosh ^, 5 = sinh ^ — 1. 
These are such that they produce the following solutions: 



(A) Plane wave: '?/'yi(x, I/, 2;) = cos(u;x) cos(a;y) cos(ci;2) 

(B) Spherically-symmetric pulse with width tq, concentrated at the origin, falling off with 
order l/r^ as r ^ cx): Vb('^) = i+r^/r^ 

(C) Toroidal solution, with radius ~ Tq and width ~ a,,, in the radial direction and cr^ in the 
vertical one: i)c{p.z) = (cosh(^) - + ^^V^"^)"^ 
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Fig. 8. Convergence of the solution error, in the I2 and loo norms, with respect to the number 
of points N, for hnear and quadratic finite elements. The plots (a), (b) and (c) correspond to the 
potentials Va, Vsir) and Vc{p,z), respectively. Prom each pair of lines, the upper one represents 
the loo norm, and the lower one the I2 norm. Each plot shows the convergence orders re, obtained by 
linear fit using the four points with highest resolution. The numbers in brackets give the convergence 
orders in the ^oo norm. 



These three potentials are chosen such that the solution for tfj is known in closed form, and 
they all differ in the way they are adapted to the underlying multi-block grid. The first 
potential, Va, represents a simple periodic wave, with the mesh not adapted to the shape of 
the wave. Second potential is supposed to model a situation where the wave is concentrated 
near the origin, in the central cubical patch of the "cubed sphere" domain. In this case, 
the grid resolution is adapted to the solution, but not the coordinate lines. Finally, for the 
potential Vc (which has toroidal shape), both the resolution (in the 6 and r directions), and 
coordinate lines (in ip direction) are adapted to the solution. 

All the test problems were solved on the same 7-patch spherical domain, with dimensions 
Rout = 10, Qc = 2.5, and fixed grid size ratios N : = I : 2. The test problems used 
the following set of parameters: for Va- uj = 0.1, for VB(r): Tq = 4, and for Vc{p, z): Tq = 8, 
o"r = 1.2, az = 4. 

It is well-known (see, for example, [56], [57], [58]) that in case of the optimal approximation, 
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the convergence rate of continuum-level error norms, \ \uh — tie| I2 and | \uh — ^^e| |oo, defined in 
a usual manner, 

\\Uh- Ue\\2= (^J^iUhix) - Ueix)ydx^ , (4.8) 
\\Uh- Ue\\oo = ^^x\Uh{x) -Ue{x)\ (4.9) 

for the standard uniform refinement, is determined by the approximation power of the finite 
element function spaces (here, Ue{x) is the exact solution, Uh{x) is its finite element approx- 
imation, and index h denotes the maximum simplex diameter in the domain triangulation 
Th). If piecewise polynomials of fixed order p are used, then the order of convergence of the 
continuum-level error norms is p + 1. 

In order to measure the error of the grid solutions, we use discrete I2 and loo norms, sampled 
at the nodes and normalized by the corresponding norm of Wg: 



maxfc \uh,k - Ue{xk)\ 
maxfc \ue{xk)\ 

The plots on figure M show the convergence of e/1,2 and e/i,oo with (A^ is inversely propor- 
tional to the maximum mesh diameter h, see section 1411!) . Note that the convergence orders 
in the I2 and l^o norms agree, which means that the pointwise convergence order is the same 
everywhere. The observed convergence order for linear elements is 2, which is supposed to 
be the case, since the order of piecewise polynomials is odd. For quadratics, we expect to 
have superconvergence, and indeed, the observed convergence rate is ~ 4. 



4-4 Testing convergence of numerical derivatives of the solution 



To set up initial data for our General Relativity evolution codes, we need not only the 
solution itself but also its first spatial derivatives, because we use a first-order formulation of 
the Einstein evolution equations. In total, our finite element solution has to be differentiated 
twice: once, when setting initial data, and one more time when computing the evolution 
equations. In this subsection we numerically study how our obtained finite element solutions 
behave under two numerical finite-difference differentiations in terms of convergence. 

If we were using completely unstructured meshes and interpolated the finite element solution 
to the multi-block grids used in our evolutions, the resulting grid solution would have an error 
of order 0{h^~^^). Two numerical differentiations in this case would take away two orders of 
convergence, leading to an unacceptable low convergence order. However, as discussed above, 
when using semi-structured grids special conditions which lead to superconvergence can be 
met, and the convergence rate of numerical derivatives improve. 
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Fig. 9. (a) Convergence of the solution error in the ^2-iiorm, for Unear and quadratic finite elements, 
for the test problem with potential Va and uj = 0.1. (b,c,d) Convergence in the /2-iiorm of the first 
(left) and second (right) numerical FD derivatives of: (b) the exact solution, restricted to the FD 
grid; (c) the numerical solution obtained using linear finite elements; (d) the numerical solution 
obtained using quadratic finite elements. In all cases several FD operators are used to compute the 
derivatives, and the resulting convergence factors (denoted by re) are shown 
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The procedure for converting the finite element solution into the grid solution is trivial: the 
value of the solution u{x) = Z^ili Ci0j(a;) at a node i is simply the corresponding nodal 
coefficient Cj. In order for the nodes to coincide with the curvilinear gridpoint in the blocks, 
we restrict our finite element solutions only to vertex nodes, and omit other types of nodes, 
such as mid-edge ones. 

In our multi-block evolutions we use new, efficient high-order finite differencing (FD) opera- 
tors satisfying the summation-by-parts (SBP) property, constructed and described in detail 
in [35]. If D is a one-dimensional differential operator, the SBP property means that for 
any grid functions u and f on a segment [a, b] with a constant grid spacing h, the following 
condition is satisfied 

{Du,v) + {u,Dv) = u{h)v{h) — u{a)v{a), 

where (■, ■) denotes the scalar product between two grid functions, defined by the SBP matrix 
S = associated with D: 

{U,V) = "^(TijUiVj 

In this paper, we use the following FD operators: -D2-1, -D4_2, D^^s, Dq^^, Dq^^ and -D8-4- 
The pair of numbers in the FD operator's subindex refiects the convergence order at interior 
points and at points at and close to the boundary. In more detail: for a generic FD operator 
Da-b, the convergence order in the interior is a and at and close to the boundaries it is 6 < a. 
The convergence of the numerical derivative in the I2 norm is at least of order 6-1-1 [59] and 
in the loo norm - it is at least of order b. 

The operators -D2-I5 -D4_27 -De-s and -D8-4 are diagonal norm (scalar product) based, as they 
satisfy SBP with respect to certain diagonal scalar product (norm) of grid functions. -D4_3 
and -Dg-s are so-called "full restricted norm" operators. They satisfy SBP with respect to 
norms which are not necessarily diagonal, but only restricted to be diagonal at the boundary 
points. The diagonal norm FD operators have several advantages compared to full restricted 
norm ones in terms of stability properties. However, they exhibit lower order of convergence 
at and close to boundary points. While the full restricted norm operators are only one order 
less convergent at the boundary, diagonal norm operators lose half the convergence order 
compared to the interior points. 

After the finite element solution is converted to a grid function and then numerically differ- 
entiated with FD, the resulting convergence factor for the numerical derivatives depends not 
only on the order of the finite element basis polynomials, but also on the order of the FD 
operator. In general one expects that FD operators of sufficiently high order will preserve the 
original convergence of the finite element solution. This is a very nontrivial mathematical 
result in super convergence theory (see Section 8.2 of [32] and references therein). In this 
section we test it for first and second FD derivatives (where by 'second FD derivative' we 
mean 'first FD derivative, applied two times', as opposed to a second-order FD derivative). 
The second FD derivative is expected be one order less convergent. 
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Table 1 

Top table: convergence orders of the solution error in the Z2-iiorm, for linear and quadratic finite 
elements, for the three test potentials. Middle and bottom tables: convergence orders of the first 
and second numerical derivatives, computed with different SBP operators. 



The results of our numerical experiments are illustrated by the figure [9] and summarized in 
tables [1] and O These tables list convergence orders in the I2 and /oo norms for each SBP 
operator, applied once and twice to three different grid functions: the exact solution, the 
numerical solution obtained with linear elements, and the numerical solution obtained with 
quadratics. Test problem is the equation 14.11 with three test potentials: Va with u = 0.1, 
Vb with ro = 4 and Vc with parameters cr^ = 10, = 5, ro = 1.5. For the potentials Va 
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Table 2 

This table displays the same type of information of the previous one, but this time in the /qo norm. 



and Vb we use a spherical domain with 7 patches, Rout = 10, Oc = 2.5, and grid size ratio 
N : iVr = 1 : 2. For Vc we use 13 patches with Rout = 20, Rmed = 7, Cc = 1.5, and grid 
size ratios N : Nr^m '■ Nr^out = 1:1:1. 

Figure [9] shows in more details some of this information, displaying log-log plots of the I2- 
norms of the errors of the solution and its first and second numerical derivatives, computed 
with various FD operators, for the problem 14. II with the potential Va- 



Several conclusions can be drawn from these tables and figure. In general, one expects the 
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convergence order of the first FD derivative to be tfie smallest between the convergence 
order of the finite element solution itself and the convergence order of the SBP operator 
used to compute the derivative. Our results support this expectation: a) the first numerical 
derivative of the exact solution converges with the order of SBP operator used to compute it. 
b) The convergence order for the first derivative of the numerical solution obtained with linear 
elements approaches 2 for all SBP operators, c) The convergence order for the first derivative 
using quadratics improves as the SBP order increases. Eventually, when the operators -D8-4 
and -De-s are used, the convergence order is either equal to the convergence of the finite 
element solution, or to the convergence order of the corresponding SBP operator. 

The numerical results also show that second numerical differentiation takes away one order 
of convergence for all three grid functions. In particular, the second numerical derivative of 
linear elements solution fails to converge in the /oo norm (see table [2]). Similarly, the second 
numerical derivative computed with the -D2-1 operator fails to converge as well. 

The most important conclusion from these results is that taking successive numerical deriva- 
tives of the grid function is only efficient for quadratics, i.e. when super convergence takes 
place. In this case, the finite element error can be smaller than finite differencing one, and 
the convergence rate of the latter is observed. For linear elements (and other elements of odd 
order), other techniques are required. One of the possible solutions is super convergent gra- 
dient recovery, when the finite element solution is differentiated and then its discontinuous 
derivative is projected back onto original finite element space (sometimes with additional 
postprocessing, see, for example, [60], [61], [50], [62] ) We do not pursue this direction, be- 
cause superconvergent gradient recovery for linear element solution will only produce no more 
than second-order initial data, while for quadratics we already have third-order convergent 
initial data without extra effort. 

4-5 Adaptive mesh refinement 

One of the main advantages of the finite element method in general and FETK in particular 
is fully adaptive mesh refinement (AMR). We have explored AMR in our semi-structured 
grids through the following strategy: we start from a multi-block triangulation, adaptively 
refine it until the error reaches some predefined level, and read off the final result at the 
original multi-block triangulation nodes [101 In doing so, we have learned that for the type 
of problems we are solving for in this paper, AMR is not necessarily the most efficient 
strategy. 

First of all, no matter how much the mesh is refined, in the end the finite element solution 
has to be restricted to the same FD grid. What matters for our purposes is the ability of 
the FD code to properly approximate the solution on this grid, and not at the refined FE 
mesh. The original FD multi-block grid is already adapted to be sufficiently fine to resolve 
all important features of the solution. 

Second, in the cases here treated both the final solution is a rather regular smooth function, 
and the domain has regular boundaries. For such cases, standard uniform refinement gives 
just as good error reduction as AMR. AMR should be advantageous, though, in cases where, 
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Fig. 10. Error in the /2-iiorm for the solutions obtained with hnear and and quadratic finite elements 
with adaptive mesh refinement. The test problem used is that one of section with potential Vq- 
The straight lines represent the fitted convergence exponents. For linear elements, at each resolution 
we started with a semi-structured multi-block triangulation and did four AMR iterations. The I2 
norms of errors for these interations are plotted with the same abscissa, corresponding to the initial 
resolution. For each resolution, the set of iterations done with quadratic finite elements is the result 
of p-refinement on the linear element meshes. 

for example, the solution has non-smooth features or the domain has non-smooth boundaries. 

Finally, the set of grid points where we sample the solution (that is, the multi-block ones) 
is very special. As already explained, the convergence order at these points is generally 
higher than one might expect. The special status of these grid points implies local symmetry 
of the finite element function spaces with respect to these points, which in turn implies 
superconvergence. AMR can easily break this symmetry and degrade both the error and the 
convergence order to the level of "ordinary" points. If we wanted to keep this symmetry and 
have superconvergence, we would have to refine at least within the entire patch, but this 
would be almost as expensive as refining the entire domain uniformly. 

Our experiments do not show particular advantage of using AMR compared to the uniform 
refinement of semi-structured multi-block triangulations. In our experiments we found that, 
in general, after a few refinement steps the AMR error saturates, while uniform refinement 
error continues to decrease. We have found that often, while decreasing the global I2 norm, 
AMR leads to an increase in the loo norm of the solution, because the local symmetry is 
broken at several grid points in the refined mesh. 

We have also tried a simple form of p-refinement. Namely, using the same set of refined 
meshes for linear and quadratic elements (see figure [TU]) . We chose several initial multi- 
block triangulation meshes with different resolutions and did four AMR iterations with 
linear elements. Then we used the same meshes to increase the order of finite elements to 
quadratics (p- refinement). It turns out that at high resolutions the meshes which give the 
best error reduction when going from linear to quadratic elements are initial triangulation 
meshes with no AMR. It happens so because the refined mesh loses the property of local 
symmetry at some of the grid points, the pointwise superconvergence for quadratic elements 
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is lost and the global L2-norm error is observed instead. 



5 Brill waves initial data and evolutions 

In General Relativity, initial data on a spatial 3D-slice has to satisfy the Hamiltonian and 
momentum constraint equations [63ll6l] . 

- K'^K.j + = (5.1) 
Vi{K'^ - g'^K) = (5.2) 

where Kij and K are the extrinsic curvature of the 3D-slice and its trace, respectively, and 
the Ricci scalar associated with the spatial metric Qij. 

Brill waves pSj constitute a simple yet rich example of initial data in numerical relativity. 
In such a case the extrinsic curvature of the slice is zero, and the above equations reduce to 
a single one, stating that the Ricci scalar has to vanish: 

^R = 0. (5.3) 

If the spatial metric is given up to one unknown function, Eq. (15. 3p in principle allows us 
to solve for such function and thus complete the construction of the initial data. The Brill 
equation is a special case of (15.31) . where the 3-metric is expressed through the conformal 
transformation Qij = ip^Qij of an unphysical metric cjij, with an unknown conformal factor 
ip. Equation (15. 3p then becomes [66] : 

(-V^ + ^R)ij = (5.4) 

where R and V? are the Ricci scalar curvature and Laplacian of the unphysical metric gij, 
respectively. 

Here we will focus on the axisymmetric case with the unphysical metric given in cylindrical 
coordinates by 

g^. = e^'^^P''^ {dp'' + dz^) + p'd^' , (5.5) 

where g(p, z) is a function satisfying the following conditions: 

(1) regularity at the axis: q{p = 0, z) = 0, |^|p=o = 0, 

(2) asymptotic flatness: g(p, z) I r^oo < 0(l/r^), where r is the spherical radius r = ^ \ 

The Hamiltonian constraint equation ( 15. 3p becomes a second order elliptic PDE, which with 
asymptotically fiat boundary conditions at r — >^ oo takes the form 
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V'^(p, z)+V{p, z)^{p,z)=Q, 



M 



V'|.-oo = l + ^ + 0(l/r^), 



(5.6) 
(5.7) 



with the potential V{p, z) given by 



We numerically solve this equation using FETK on the 13-patch multi-block spherical 
domain described in section 14.11 (see figure [7]). We use domain parameters Rout = 30, 
R-med = 7, ac = 1.5, and grid dimension ratios : N^^inner '■ Nr^outer = 2 : 3 : 12. Our 
low- medium-high resolution triple is = 32, = 36 and = 40, except for pointwise 
convergence tests on the x-axis (see figure [T3l) . where we use A^ = 16, A^ = 24 and A^ = 36 
(since they all differ by powers of 1.5). 



We impose Robin boundary conditions, as in equation (14.51) . The weak form (14.61) and bilinear 
linearization form (14. 7p for this problem are given in section 14. 2[ Since first order elements 
lead to unacceptably low convergence orders for most general relativistic applications, from 
hereon we restrict ourselves to quadratic ones (which should give fourth order convergence 
if superconvergence is exploited) . 




(a) 



(fa) 



Fig. 11. Potentials for the two types of Brill waves considered: Holz' (a) and toroidal (b) forms. 

We work with two specific choices for g(p, z): 

(a) Holz' form [67]: qnip, z) = anp'^e'^^ , with amplitude an = 0.5; 

(b) toroidal form: qt{p, z) = atp^ exp ( —is^:^ with amplitude at = 0.05, radius po = 5, 



width in p-direction = 3.0 and width in 2;-direction = 2.5. 

Holz potential is chosen for historical reasons. It is suitable for evolutions using the cubed 
sphere system, because initially the wave is concentrated near the origin (at the central 
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patch), and later decays into a sequence of spherical waves. The cubed sphere grid resolves 
the wave well both at the initial moment, and during the subsequent evolution. For the 
toroidal potential, the cubed sphere system is adapted even better, since it efficiently removes 
the dependence on one of the angular coordinates, (/?. 

5.1 Importing the initial data into QUILT 

Once we have constructed the initial data sets, we analyze and evolve them by importing 
them into QUILT. For evolutions we use the Generalized Harmonic (GH) first order sym- 
metric hyperbolic formulation of the Einstein's equations introduced in [68], which features 
exponential suppression of short- wavelength constraint violations; our multi-block implemen- 
tation is described in [37] . 

The set of evolved variables in this system includes the 4-metric Qab, its spatial derivatives 
^iab = diQab and quantities Hat = —f^dcgab, where is the unit normal vector to the spatial 
slice. 

To set up the initial data, we first compute the 3 + 1 quantities and then convert them to 
GH variables. The 3-metric gij{x,y,z) is computed from the conformal factor 2), in 

Cartesian coordinates, using the expressions (which follow from 15.51) : 



where p = a/x^ + y'^. Then we construct the rest of the evolved variables, including the gauge 
source functions Ha = -g^Ta,bc (here Ta^bc = \{,dbgac + dcgab - dagbc) are the Christoffel 
symbols). The extrinsic curvature is assumed to be zero, Kij = 0, we also use unit lapse 
a = 1, zero shift (3i = 0, and zero time derivative of lapse and shift: dta = dtPi = 0. 

5.2 Convergence of initial data and Hamiltonian constraint 

To estimate the quality of our initial data, we evaluate the Hamiltonian constraint violation 
using the SEP operators, which is in fact equivalent to an independent residual evaluation 
for the Brill equation (15. 6p . As is the case with just numerical derivatives, we find that the 
magnitude and convergence order of the Hamiltonian constraint violation depends on both 
the order of the SEP operator, and the order of finite elements. In total, computing the 
Hamiltonian constraint involves two numerical differentiations, therefore it has to converge 
with the same order as second numerical derivatives. We can see this convergence rate for dif- 
ferent SEP operators in figure [T2l Table [3] summarizes those results, along with the expected 
convergence orders for second numerical derivatives. 

Consistent with what one expects, for finite difference operators of sufficiently high order 
the constraints converge with the same order as the finite element solution itself, which 
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should be 3 + (T for some < cr < 1 (depending on the level of superconvergence obtained). 
Similarly, below in section 15.31 we will show that the extracted gravitational waves have a 
similar convergence order. 
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Fig. 12. Convergence of the Hamiltonian constraint in the Z2-norm for Brill waves of Holz (left) 
and toroidal (right) type, computed using quadratic finite elements and different numerical SEP 
operators. 
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Fig. 13. (a) 1-d cut through the Brill wave conformal factor V'/me for problem (a) along the x-axis. 
(b) Errors ipcoarse - ipmedium, tpmedium -ipfine, and poiutwisc Convergence order on the 1-d cut along 
the X-axis. 

It is not difficult to see why the convergence order is lower for problem (a). This is the order 
we might expect for quadratic elements on completely unstructured meshes, without super- 
convergence. The reason that no superconvergence is observed in this case is the following. 
We have chosen the domain and the width of the Brill wave in such a way that around the 
boundary of the inner cubical patch the solution varies significantly (see figure [T^ . Recall 



24 



that the size of that patch is ac = 1.5 and the width of the gaussian in the function q is 1. 
But this is exactly the place where the local symmetry property is violated most (especially 
at the corners of the cube), and conditions for superconvergence are the least favorable. Ev- 
erywhere else the solution varies very slowly and the error is small compared to the error at 
the boundary of the central cubical patch. With increasing resolution this error dominates 
in both I2 and /qo norms. The situation is better for the problem with toroidal potential, 
because most of the variation of the solution is located outside the central cube (recall that 
the radius of the toroidal wave we use is tq = 5); though this is still in the inner-patches 
region (the radius of the spherical boundary between inner and outer patches is r^ed = 7). 

5.3 Multi-block evolutions 

We now demonstrate that our approach for generating initial data on multi-block grids using 
finite element methods can be successfully used in practice in fully nonlinear relativistic 
simulations. We do so by presenting results of multi-block evolutions of the Holz set of Brill 
initial waves constructed above. 

In the notation of [HS], we fix the damping parameters of the GH formulation to 7o = 72 = 
1. We used the SBP operator -De-s for spatial differentiations, a 4-th order Runge-Kutta 
time integrator with adaptive time stepping, and maximally dissipative outer boundary 
conditions. 

The Brill wave amplitude an = 0.5 is in the subcritical regime (the critical value is around 
ttcr ~ 4.85 [SZl). As a result the wave, initially concentrated near the center, dissipates and 
leaves the domain after a while. Figure [T3] shows a convergence plot in time for the Hamil- 
tonian constraint during such evolution. We see that the Hamiltonian constraint converges 
with a factor of 2 — 3 in the I2 norm. This has to be one order less convergent than the 
solution itself, therefore we anticipate the solution to converge with a factor of 3 — 4, which 
is in agreement with the pointwise convergence order of the conformal factor at the initial 
time (figure [T^ . 

We compute gravitational waveforms using a generalized Regge-Wheeler-Zerilli formalism, 
along the lines of [37] and extending that reference to the even parity sector (details about 
that extension will be presented elsewhere). The four- dimensional spacetime metric is de- 
composed into a spherically symmetric "background" plus a small perturbation [69|70] . The 
background part of the metric is identified with the Schwarzschild solution and the radiation 
content with the difference between the numerically computed solution and the background. 

The background metric is written as 

ds"^ = gp^{t, r)dxPdx^ + f{t, r)gABdx^dx^ . (5.8) 

with the four-dimensional background manifold split as the product space of a two-dimensional 
one }A endowed with coordinates x'' (with p, r = 0,1 usually denoting the time and radial 
coordinate) and a unit 2-sphere S*^ with coordinates x^ (with A = 2, 3 commonly taken as 
the 9 and polar spherical coordinates). Here gp^- is the metric of the manifold }A and is 
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1 2 3 4 5 10 15 20 25 

Fig. 14. Convergence of the Hamiltonian constraint as a function of time during a multi-block 
evolution of the Brill wave initial data generated with finite element methods. 

a positive function of x^. If using the areal radius as a coordinate we have f = r. The metric 
of the 2-sphere is taken to be qab = diag(l, sin^ 9) in polar spherical coordinates. 

The metric perturbation is decomposed in terms of scalar, vector and tensor spherical har- 
monics [7T]l72|37j . The decomposition naturally splits the different (£, m) modes into even 
(— 1)^ and odd (—1)^"''^ parity under reflections about the origin. The two different parities 
are handled separately. Odd and even-parity perturbations are described by the Regge- 
Wheeler [69] and Zerilli [70] functions, respectively. 

The dominant modes in the evolutions of the Brill data constructed above are the £ = 2,4 
even parity, axisymmetric ones (see figure fT5l) . Figures fT6] and [T71 display the corresponding 
Zerilli functions and their convergence behavior, extracted at a radius r^. = 12.75. The 
observed convergence factors for the i = 2 and i = A modes are around four and three, 
respectively, which are consistent with the convergence factors from quadratic elements with 
superconvergence for the initial data and the -Dg^s SBP operator and fourth-order Runge- 
Kutta for the evolution. 
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Fig. 15. Nonzero components of gravitational radiation i = 2 and £ = 4, m = 0, extracted at radius 
re = 12.75. 
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Fig. 16. Zerilli function (left) and self-differences (right) for the {i 
according to fourth order convergence. 
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Fig. 17. Zerilli function (left) and self-differences (right) for the {£ 
according to third order convergence. 



4, m = 0) mode, scaled 



Final remarks 



In this paper we followed a finite element approach for generating initial data satisfying the 
Einstein constraint equations on semi-structured, multi-block three-dimensional domains. In 
section H] we used semi-structured multi-block triangulations to solve for some test prob- 
lems with known closed-form solutions. The obtained linear and quadratic finite element 
solutions were then restricted to the multi-block grid, and their convergence, as well as the 
convergence of their first and second derivatives, was evaluated numerically using indepen- 
dent high-order finite difference operators satisfying summation by parts (SBP). While the 
linear elements solution showed usual 2-nd order convergence (unacceptably low for many 
relativistic applications), for quadratic elements we obtained superconvergence with order 
3 + cr (with < cr < 1) on the grid, due to the approximate local symmetry at the mesh 
vertices with respect to the local inversion of the multi-block triangulations. 

Initial data for a first-order formulation of the Einstein equations involves first derivatives 
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of the solution of the constraints equation. Computing the constraints or right-hand sides 
of the evolution equations requires taking derivatives twice. In subsection 14.41 we analyzed 
convergence of the first and second numerical derivatives, taken with different high-order SBP 
operators. For quadratic elements, the first numerical derivative was observed to converge 
with either the super convergence order 3 + a, or the order of SBP operator. The latter is 
a transient error behavior and happens when the FE error is still smaller than the error of 
numerical differentiation. 

In subsection 14.51 we discussed three factors which make adaptive mesh refinement (AMR) 
unnecessary and/or less efficient for the problems here considered when compared to global 
refinement: the fact that the multi-block grid is already tailored to resolve fine features of the 
solution, the need to restrict the finite element solution to the same grid, and the superconver- 
gence properties of the quadratic elements solution. Because we lose superconvergence when 
using completely unstructured meshes, adaptively refining the solution sometimes makes the 
errors larger (see figure [TD] for an example). However, we also noted that AMR would likely 
become advantageous for other problems with more singular solution features. 

Finally, in section [5] we presented numerical experiments with Brill waves. The constraint 
equations in this case reduce to a single elliptic one (15. 6p on the conformal factor ip, which 
has to be differentiated once to obtain the full set of initial data variables in the generalized 
harmonic formulation (subsection EH]) . Subsection 15. 2l presented a convergence analysis of the 
initial data and Hamiltonian constraint, and confirmed that the initial data computed with 
quadratic finite elements shows the desired order of convergence > 3 (see table [3]). Finally, 
in section 15.31 we demonstrated stable, > 3-rd order convergent multi-block evolutions of 
subcritical Brill waves with finite differences, summation by parts operators, and extracted 
the first two dominant radiation modes from the numerical solution. 

This paper shows that generating initial data on semi-structured multi-block triangulations 
using finite element methods is a feasible approach which works well in practice. Future work 
might include adding higher order and/or spectral elements to this approach. 
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